Coagulation reactions in low dimensions: Revisiting subdiffusive A + A reactions in 

one dimension 
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We present a theory for the coagulation reaction A + A —* A for particles moving subdiffusively 
in one dimension. Our theory is tested against numerical simulations of the concentration of A 
particles as a function of time ( "anomalous kinetics" ) and of the interparticle distribution function 
as a function of interparticle distance and time. We find that the theory captures the correct 
behavior asymptotically and also at early times, and that it does so whether the particles are nearly 
diffusive or very subdiffusive. We find that, as in the normal diffusion problem, an interparticle gap 
responsible for the anomalous kinetics develops and grows with time. This corrects an earlier claim 
to the contrary on our part. 
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I. INTRODUCTION 



The traditional laws of mass action that describe the 
time evolution of the macroscopic global concentrations 
of reactants and products in chemical reactions assume 
that the system is well stirred and therefore spatially ho- 
mogeneous. However, there are many situations when a 
reactive system is not well mixed; in that case one must 
deal with local concentrations and account for the ef- 
fects of spatial inhomogeneities on local reaction rates. 
Spatial variations in concentrations of reactants lead to 
changes (often called "anomalies") in the time depen- 
dences of the spatially averaged macroscopic concentra- 
tions. One often encounters this situation when diffu- 
sion is the only mixing mechanism, particularly when 
diffusion is the rate limiting step for a reaction to oc- 
cur. The inefficiency of diffusion as a mixing mechanism 
becomes more pronounced with decreasing dimensional- 
ity of the system, and it is therefore commonly accepted 
that diffusion-limited reactions in constrained geometries 
exhibit kinetic "anomalies" [lj. An exact description of 
diffusion-limited reactions in the face of nonuniform spa- 
tial distributions of reactants typically requires an infinite 
hierarchy of correlation functions to properly incorporate 
spatial correlations. In practice, such hierarchies are of- 
ten truncated at the first or second level, giving rise to 
well-known reaction-diffusion equations whose solutions 
are sometimes fairly satisfactory (and sometimes not) in 
capturing the principal deviations form the laws of mass 
action 0,0]. 

The judgment as to the success or failure of approxi- 
mate reaction-diffusion models has mostly relied on com- 
parisons with numerical simulations. There has always 
been a quest for exact analytic solutions against which 
approximate reaction-diffusion theories could be tested, 
but there have been very few successes, among them the 
coagulation (A + A — ► A and A + A ;=± A) and annihila- 
tion (A+A — > 0) reactions . Exact solutions in these 
cases have been possible because, if instead of focusing on 



the concentration of reactants, one focuses on the evolu- 
tion of empty intervals (coagulation reactions) or on the 
number parity (even or odd) of particles in an interval 
(annihilation reactions), one arrives at exactly linear dif- 
fusion equations. These solutions have provided a wealth 
of information against which to measure the approximate 
solutions obtained from more standard approaches for 
these particular reactions. Unfortunately, the interval 
approaches are not generalizable to other reactions. 

In the past few years there have also been attempts to 
understand chemical reactions in low dimensions when 
the reactants move subdiffusively. This has been a par- 
ticularly interesting subject in view of the many systems, 
mainly biological, in which such reactions occur in a com- 
plex or constrained environment that does not even per- 
mit ordinary diffusive motion of chemical species. All 
the difficulties encountered in diffusive systems are ex- 
acerbated in this case because motions are even slower, 
a fact that turns out to have profound consequences on 
spatial as well as temporal correlations Q and to cause 
much greater difficulties in both numerical simulations 
and analytic attempts. One approach to this problem has 
been to essentially adapt the existing reaction-diffusion 
models by modifying the diffusive description to a subdif- 
fusive one involving fractional diffusion operators. These 
descriptions are phenomenological in nature, and much 
work remains to be done to understand how one might ar- 
rive at a mesoscopic description from microscopic consid- 
erations. In particular, there are many questions about 
how to describe the reaction terms in subdiffusive sys- 
tems, and there is clear evidence that the problem is 
usually not even separable into a simple sum of a term 
describing the motion of the particles and one describing 
the reaction. 

The subdiffusive coagulation and annihilation prob- 
lems seemed to offer a parallel opportunity for exact so- 
lution if one again concentrated on the properties of in- 
tervals. We followed this route in our earlier work [f|, 
and formulated what we thought to be exact subdiffu- 
sive equations for the intervals from whose solution one 
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could calculate the particle concentrations and interpar- 
ticle separations, as had been done in the diffusive prob- 
lem. These solutions led us to the conclusion that in- 
terparticle gaps do not occur in the subiffusive problem, 
a result that seemed reasonable for particles that move 
sufficiently slowly. However, subsequent numerical simu- 
lations indicated that a gap does develop no matter how 
subdiffusive the particles, and this led to a reassessment 
of the assumptions of our original theory. In this pa- 
per we present this numerical evidence, a description of 
the difficulty with the original theory, and a new the- 
ory which, albeit still approximate, seems to capture the 
correct behavior to a very high degree of accuracy. Also, 
we recently developed a mean field theory for this prob- 
lem Q but it is only valid for dimension 3 and higher. 
We concentrate on the coagulation problem, although a 
similar approach may be helpful for the annihilation re- 
action. 

In Sec. |TT] we describe the difficulties with our previ- 
ous ID theory |5|. Section |TTT] presents our new theory. 
Section ITVl discusses some special interesting issues that 
arise in this problem, and in Sec. [V] we present the numer- 
ical evidence to support our theory. We conclude with a 
summary in Sec. IVII 



II. TRADITIONAL THEORY 

The interparticle distribution function method focuses 
on the "empty interval" function E(x,t), defined as the 
probability that an interval of length x is empty of parti- 
cles at time t. From E(x, t) one obtains the concentration 
of particles, 



c(t) = - 



dE(x,t) 



Al- 



and the interparticle distribution function, 



p(x,t) 



1 d 2 E(x,t) 
c(t) dx 2 



(1) 



(2) 



which is the probability density that the first A to be 
found to one side of a given A at time t is a distance x 
away. The question of interest is how to determine the 
empty interval function. 

If the particles undergo normal diffusion, the probabil- 
ity density for finding a particle A at y at time t in the 
absence of reactions obeys the diffusion equation 



0_ 

Of 



P(v,t) 



(3) 



In the presence of the coagulation reaction, the empty 
interval function focuses on the diffusive motion of the 
particles at each end of an empty interval and one readily 
arrives at the exact equation 



d_ 

Of 



E(x,t) 



2D^E(x,t). 



(4) 



The equation is easily understood from the fact that the 
empty interval dynamics is the same as that of the indi- 
vidual particles but with double the diffusion coefficient 
to reflect the relative motion of two diffusive particles. 
This readily tractable equation, together with appro- 
priate boundary conditions, exactly solves the diffusive 
A + A — ► A problem in one dimension. 

A standard approach for the description of subiffu- 
sive processes starts from the continuous time random 
walk (CTRW) formalism, in which a walker jumps from 
one site on a lattice to another in consecutive steps as 
time proceeds 0, II] . Both the jump distances n and 
times t between jumps are random variables drawn from 
a probability distribution function <j)(n,t). If the jump 
distances and jump times are independent random vari- 
ables, then this distribution function is simply a product, 
<p(n,t) = w(n)ip(t). "Normal" CTRWs are associated 
with distributions w(n) and ip(t) that have finite first mo- 
ments, and the scaling limit that leads from the random 
walk to the diffusion equation is well known. One way 
to obtain a subdiffusive process is for the waiting time 
distribution ijj(t) to be heavy-tailed, i.e., ip{t) ~ t -1-7 
with < 7 < 1 for long times, so that the mean waiting 
time between jumps diverges. In this case a number of 
scaling approaches can be found in the literature, with 
a particularly helpful discussion in [§] . In the absence of 
reactions, in the continuum limit with a particular scal- 
ing one arrives at a "fractional diffusion equation" for the 
evolution of the probability density P(y, t) of a subdiffu- 
sive particle, 



fl B 2 
— p(„ A - ^n 1 -"" 

dt 



P(y,t)= oD r a K a —P(y,t), 



where qD\ a is the Riemann-Liouville operator, 
° Dt P{y ' t) = T{a)d-tJ dT (t-r)i-« ' 



(5) 



(6) 



and K a is the generalized diffusion coefficient. The mean 
square displacement of the A particle for large t that 
follows from this evolution equation is 



2K n 



r(l + a) 



r 



(7) 



which reduces to the ordinary diffusion result when a — 1 
(Ki = D). In our earlier work we argued that the same 
reasoning that led from Eq. ([3]) to Eq. would also lead 
from Eq. (O to the interval equation 

g- t E(x,t)= vD]- a 2K a — E(x,t), (8) 

and then calculated the particle concentration and in- 
terparticle distribution function from the solution of this 
equation with appropriate boundary conditions. 

The difficulty with this reasoning is the fact that in the 
subdiffusive problem we must face the issue of aging [To| , 
which we have not done above. To describe the issue, we 
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again turn to a CTRW point of view of the problem. 
Consider first the diffusion problem, where the contin- 
uum limit leads to a diffusion equation. In this limit, in 
which the mean time between steps and the step length 
go to zero in an appropriate way, one arrives at the dif- 
fusion equation. In particular, in the diffusive problem 
the steps are sufficiently frequent that one need not keep 
track of extremely long sojourns of a particle at any one 
site, and what transpired in the past is quickly forgotten. 
A similar formulation of the subdiffusive problem involves 
waiting time distributions with an infinite mean time be- 
tween steps. Even worse, because of aging, the waiting 
time distribution for a particle to take its next step at 
time t having arrived at its current location at time t' is 
no longer simply a function of the difference t — t' but of 
each time separately. To understand the effects of these 
complications, suppose that an observer looks at the sys- 
tem at time t and sees an empty interval of size x at that 
instant. For normally diffusive particles, the evolution of 
the size of this interval does not depend on the time at 
which the interval was first created. However, in order 
to predict the evolution of this gap for subdiffusive par- 
ticles the observer must know how long each of the two 
particles at the ends of this interval have been at that 
location. Due to aging, the evolution depends on the 
times t — ti and t — t r , where ti and t r are the times at 
which the left and right particles jumped to the locations 
seen by the observer at time t. Moreover, even if the two 
particles arrived at this location at the same time, i.e., 
ti = t ri the shortening and lengthening of this gap is not 
correctly described by Eq. ([5]),as we will show in Sec. IIIII 
[cf. Eqs. ©-(El)]. 

III. THEORY REVISITED 

To find a more appropriate description for the evolu- 
tion of the empty interval, we introduce two hypotheses. 
As we state each hypothesis, we provide an argument as 
to its approximate nature. 

We start by forgetting the reaction for a moment, and 
simply consider the motion of (mutually transparent) A 
particles on an infinite line. We define V a (x,t) as the 
probability density for x to be the distance between two 
subdiffusive particles at time t with the initial condition 
V a {x,0) = S(x), i.e., 

dyP a (x-y,t)P a (y,t). (9) 

-oo 

Here P a {y, t) is the propagator of the subdiffusion equa- 
tion, i.e., the solution of Eq. ((5|) with the initial condition 
P a (y,0) = S(y). In Fourier space the solution in closed 
form is the Mittag-Leffler function, 

P a {q,t)=E a (-K a q 2 t a ). (10) 

The Fourier transform of Eq. ^ then tells us that 

V a (q,t) = (P a (g,t)) 2 = [E a (-K a q 2 t a )] 2 , (11) 



and correspondingly, 




dqe-^ [E a {~K a qH a )} 2 



1 f°° 

= - dqcos(qy)[E a (-K a q 2 t a )] 2 .(12) 

T Jo 

We do not know the evolution equation for V a (x,t), 
i.e., we do not know the operator 5" Q such that 
■SaPa{x, t) = 0, except for a = 1. In this case the Mittag- 
Leffler function reduces to an exponential, E\(— z) = 
exp(— z), so that Vi(q,t) = exp(— 2K\q 2 €) and the corre- 
sponding function Vi(x, t) is the solution of the diffusion 
equation ^ but with double the diffusion coefficient, 

d d 2 

—T 1 (x,t) = 2D—T 1 (x,t), (13) 

with K\ = D. However, for a < 1 the Fourier transform 
of the solution of the subdiffusion equation with double 
the subdiffusion coefficient, 

Q- t Q<*(x,t) = Q D\- a 2K a ^Q a {x,t), (14) 

is the Mittag-Leffler function of twice the argument in 
Eq. (HDJ), 

Q a (q,t) = E a (~2K a q 2 t a ). (15) 

Clearly, except for a = 1, 

[E a (-K a qH a )] 2 ? E a (~2K a q 2 t a ), (16) 

so that V a 7^ Qa ) and the operator $ a is not straightfor- 
wardly obtained from the subdiffusion equation, 

9a * | - oD\~ a 2K a ^. (17) 



A. Hypothesis 1 

The central hypothesis of our new theory is that 
the (unknown) equation that describes the evolution of 
V a (x,t) is the same as the equation that describes the 
evolution of the empty interval function E(x,t), i.e., that 

3 a E(x, t) = 0, 0<x<oo (18) 

This hypothesis is in general an approximation because 
the distribution of particles is not expected to be the 
same in the presence and absence of the reaction (al- 
though in higher dimensions it is 6]). While there is no 
memory of prior reaction events in the diffusive prob- 
lem, in the subdiffusive case this memory persists and 
may lead to a distortion of the distribution relative to 
that of the particles when there is no reaction. We only 
know how to assess the severity of this approximation via 
comparison of results with numerical simulations Q (see 
Sec. E| . 
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How does this hypothesis help us find the empty inter- 
val distribution when in fact we do not know the operator 
3a? It helps us because we know the solution V a (x,t), 
which is then the Green function or propagator for the 
empty interval distribution. In other words, while we do 
not know the equation of evolution for E(x,t), we have 
sufficient information to construct the function E{x, t) 
itself explicitly. The interval function is subject to the 
additional conditions 



E(0,t) = 1, 
E(oo,t) = 
E(x,0)=f(x). 



(19) 



The first boundary condition is the probability that an 
interval of vanishing width is empty. As in the diffu- 
sive problem, this probability must be unity. The second 
simply recognizes the existence of particles at any time 
t, and the third is determined by the initial distribution 
of empty intervals. 

In order to construct the solution E(x,t) from our 
knowledge of the propagator V a (x,t), we introduce our 
second hypothesis. 



B. Hypothesis 2 

The second hypothesis of our theory is that the oper- 
ator 3a is linear. We know this operator to be linear 
when a = 1 , but it is most likely an approximation when 
a < 1, although our lack of information about 3a makes 
it difficult to assess. While evolution equations describ- 
ing various quantities in some other reaction-subdiffusion 
problems are in fact not linear, the connection between 
those systems and the one considered here is not clear. 

In any case, under this hypothesis we can split E(x, t) 
into two parts, each of which satisfies conditions whose 
symmetry properties allow for the advantageous use of 
our knowledge of the propagator. In particular, we write 



E(x, t) = E T (x, t) + E A (x, t), 
where the individual pieces satisfy the following: 
3 a E T (x,t) = 0, 0<a;<oo 
E T (0,t) = 0, 
E T {oo,t) = 
E T (x,0) = f(x) 

and 

3 a E A (x, t) = 0, 0<x<oo 
E A (0,t) = l, 
E A (oo,t) = 
E A (x,0) = 0. 



(20) 



(21) 



(22) 



The superscripts T and A denote "transient" and 
"asymptotic," respectively, for reasons that become evi- 
dent below. 



C. Asymptotic (Long Time) Results 

In general, the solution E(x,t) and its derivative ob- 
servables depend on the initial distribution fix). How- 
ever, since the evolution of the interval function is in some 
sense necessarily subdiffusion-like, we expect that the 
given boundary conditions lead to a decay of E (x,t), 
i.e., E T (x,t) — > as t — ► oo, while E A (x,t) can not de- 
cay. Thus at sufficiently long times the behavior of the 
interval function will be dominated by that of E A . One 
must keep in mind that the decay of E T may be slow 
(especially when comparing with numerical simulations). 
In particular, whereas diffusive modes decay exponen- 
tially, subdiffusive modes typically decay only as t~ a for 
a < 1. Nevertheless, this provides a helpful clement 
if one is interested in the asymptotic behavior because 
E A {x, t) does not depend on the initial distribution and 
can therefore be pursued once and for all. We therefore 
focus on it first. 

The solution for E A (x,t) can be obtained by the 
method of images. For this purpose, we consider the 
related problem 



$ a £ A {x,t) = 0, 
£ A {-oo,t) = 2, 
^(oo,*) = 
£ A (x,0)=2- 



-oo < x < oo 



(23) 



26(x), 



where 9{x) is the Heaviside step function. The symmetry 
of the problem immediately leads to the conclusion that 
£ A (0,i) = 1. Therefore the solution of this problem for 
x > is just the solution E A (x,t) of Eq. (|2"2"|) , that is, 
£ A (x,t) — E A (x,t) for x > 0. On the other hand, we 
can write 

r° 

£ A (x, 0) = 2 - 29{x) = 2 1 dyS{x - y). (24) 

J — oo 

Since we know that the Green function of 3a is V a {x, £), 
and since we assume that 3a is linear, the solution of 
Eq. is 



£ A (x 1 t) = 2j ,UiP,A.r -u.l) Vir>) 
or, equivalently, 

£ A {x,l) = 2\ ,///P„ (//•/)• m\ 

Therefore 

E A (x,t) = 2l dyV a (y,t) for x > 0, (27) 



with V a (y,t) given in Eq. (TT21 . It is noteworthy that 
E A (x, t) depends on x and t only via the variable 



x, 



(28) 
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where 



dy [E a (-y 2 



(29) 



The changes of variables w = \JK a t a q and u 
(ca/VK a t a )y immediately lead to 

E A (x,t) = E(z) 
2 



du / dw cos ( ^— ] [E a (— w 2 )] 2 . 

(30) 

The particle density as a function of time is obtained 
from the interval distribution function via Eq. |T]). The 
asymptotic density lim^oo c(t) then is 

dE A (x,t) 



c(t) = 



dx 

X 



2V a (0,t) 



x=0 



(31) 



When a = 1 this reduces to the familiar result c(t) — 
[2^Dt)~ 1 / 2 . 

At these long times, using the asymptotic portion 
E A (x,t) of the solution E(x,t) in Eq. © we obtain 



c(t)p(x,t) 



d 2 E A (x,t) = _ 2 dV a (x,t) 



dx 2 



dx 



2 f°° 

- dqq sin(qx) [E a {-K a qH a )\ 



2 r , . / yx 

/ dy y sm ■ 



[E a (-y 2 )y. 

(32) 



Dchning p(z)dz — p(x,t)dx where z is defined as above, 
we find that for long times 



Kc z a Jo \c a J L 



(33) 



D. Results Valid for All Times for Initial Poisson 
Distribution 

To find the particle density and interparticle distri- 
bution function for all time requires the solution of the 
system (|2ip . which in turn requires specification of an 
initial distribution. We choose an initial Poisson distri- 
bution for this analysis, 

E r {x, 0) = f{x) = E(x, 0) = e- C0 *, (34) 

where Co = c(0) is the initial concentration of particles. 
As before, we formulate a related problem, 



5 a £ T (x,t) =0, 
£ T (-oo,t) = 0, 
£ T {oo,t) = 

£ T (x,0)=f(x) 



-OO < X < 00 



(35) 



By symmetry we deduce that £ T (0,t) = 0, so that the 
solution of this problem for x > is just the solution 
E T {x,t) of Eq. d35j), that is, £ T {x,t) = E T (x,t) for x > 
0. We can write 



/OO 
dy£(y,0)S(x-y). 
-OO 



(36) 



Then our knowledge of the Green function for $ a and 
our hypothesis that this operator is linear immediately 
leads to the solution 



x > 
x < 0. 



£ T {x,t) 



dy£(y,0)V a (x-y) 

-OO 

f dye^Vaix - y) 

J — OO 

dye- c °yV a (x-y) 



(37) 



E 1 (x,t) 



for x > 0. 



With a bit of manipulation and explicit insertion of 
Eq. (fl"2"J) we finally obtain for x > 0, 



(38) 



E {x, t) = — 2 sinh(c x) / dye- C0V V a (y) 

J X 

+ 2e- cox f dycosh(c y)V a (y). 
Jo 



The full interval distribution is then the sum of (j2~7| 
(valid for any initial condition and determinative of the 
asymptotic behavior) and (|3"8|) (explicitly calculated here 
for a Poisson initial distribution and going to zero asymp- 
totically). 

As before, the particle density as a function of time 
is obtained from the interval distribution function via 
Eq. (!]). From Eq. we calculate 



dE T (x,t) 



p OO 

= -2V a (0,t)+2c dye- c °yV a (y,t) 7 

0X x=0 JO 

(39) 

from which upon addition of Eqs. (|31[) and (|B"9l and use 
of Eq. (TT3]) it follows that 



c(t) 



dy- 



1 + (y 2 /c 2 K a t<*) 



(40) 

Note that this result, valid for all times, reduces to (j3"T|) 
when t — > oo . 

The interparticle distribution function follows from 
Eq. (12]). To add to our previous asymptotic result, we 
note that 



8 2 E J (x,t) 
dx 2 



3V 

2^-2 C gsinh( Co x) / (///< "•'"•'T., (;/./! 

+2cle~ c ° x 



dycosh.(coy)T> a (y,t). (41) 
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Introducing Eq. (fT2|) into (|4Tj) . one finds, after some ma- 
nipulations, 



8 2 E T (x,t) =2 dV a , 2 



9x 2 



9a; Trir^i J I + y 2 /c 2 K a t a 

2 



Upon addition of this contribution and the asymptotic 
one obtained earlier, we finally have 



c(t)p(x, t) 



dy 



wK a t a J 1 + y 2 /c 2 K a t<* 



(43) 



This expression reduces to Eq. (|32p when £ — > oo. Only 
in the asymptotic limit is it possible to write p(x,t)/c(t) 
as a function of the single combined variable z = c(t)x. 



IV. A CONUNDRUM AND SOME CHOICES 

In the next section we compare our theory to numeri- 
cal simulation results and, as we shall see, the comparison 
is on the whole very successful. And yet we know that 
the theory is approximate - indeed, we introduced two 
hypotheses that are surely not exact. We note here an 
additional related conundrum which exhibits itself (al- 
beit only weakly even when a is small) in the numerical 
simulations of the A + A — ► A reaction (but not in the 
A + A — > problem) . In the discrete version of the prob- 
lem used for the simulations, a reaction A+A — ► A occurs 
when an A particle steps onto a site already occupied by 
another A particle. One then has to decide which of the 
two particles is the one that is removed from the system, 
the one that was there ("kill" rule) or the one that just 
stepped onto the site ("no-kill" rule). The choice could 
vary with each reaction event. In the diffusive problem 
(a = 1) the choice does not matter. In a subdiffusive sit- 
uation, however, the choice does matter since the prop- 
erties of the subsequent random walk of the survivor, 
including the probability that the walker continues to re- 
main at that site, depend on its age. In particular, if the 
survivor is the one that arrived at the site first then the 
probability that it will remain at that site is greater than 
if the survivor is the new arrival. The implication is that 
the mesoscopic results such as the interparticle distribu- 
tion function depend on the microscopic reaction rule. 
On the other hand, we will show that the differences in 
the observable quantities are small even for a consider- 
ably smaller than 1, but the agreement with simulations 
is a bit better using the "kill" rule. This is as expected 
since this rule resets, so to speak, the "clock" of the sur- 
viving particle following a reaction event and is thus in 
some sense closer to the assumptions inherent in Hypoth- 
esis 1. In most of our simulations we use the "kill" rule. 



The second choice we must make in our simulations 
is the form of the distribution of waiting times between 
steps. In most cases we use a Pareto-type waiting time 
distribution, 



M*) = 



*o(l + t/to) 



1+a ' 



(44) 



(42) while in some we use a Mittag-Leffler-based distribution, 



lpML(t) = --E a [-(t/t ) a ] 



(45) 



(in our simulations we set to = 1). They can both be 
used to describe the asymptotic behavior of subdiffusive 
random walkers. While either form is mathematically 
acceptable, jWl is preferable at short times specially 
as a — > 1 [11]. In most of our simulations, times are 
long and a is not close to unity, so the choice is not a 
central issue. Nevertheless, we address this issue here 
even before presenting our simulation methodology and 
results because while both of these waiting times do lead 
to a fractional diffusion equation at long times, the choice 
of the waiting time distribution determines the value of 
the generalized diffusion coefficient K a . For the Pareto 
case K a p = l/2r(l — a), while the Mittag-Leffier form 
(in our simulation units) leads to K u ,ml = 1/2. In the 
discussion of our results this is the only difference be- 
tween theoretical results labeled by P and those labeled 
by ML. The Pareto generalized diffusion coefficient di- 
verges as a. — * 1. This behavior is symptomatic of other 
problems associated with this waiting time distribution 
in this limit for the calculation of quantities that explic- 
itly depend on K a , and is a strong motivator for the 
introduction of the Mittag-Leffier form. We will only use 
the Pareto distribution in unproblcmatic regimes of a, 
where these difficulties are not an issue. 



V. NUMERICAL EVIDENCE 

In this section we present numerical simulation re- 
sults to test the adequacy of our theory. Specifically, 
we present simulation results for the particle density and 
for the interparticle distribution function. We begin by 
briefly describing our numerical simulation methodology. 

We proceed via the following steps in our numerical 
algorithm. First we generate a Id lattice. We then gen- 
erate an escape time for each particle chosen from the 
given waiting time distribution. We choose the particle 
with the smallest waiting time. This particle jumps to 
one of its two nearest neighbors. If the destination site 
is empty, we update the waiting time for the arriving 
particle. If the destination site is occupied, the particles 
coagulate into a single particle. This is where we have 
to specify whether the event is of the "kill" variety or of 
the "no-kill" variety, as described earlier. Next, we look 
again for the particle with the smallest waiting time. The 
other particles that have not yet moved simply continue 
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evolving according to their own internal clocks. The par- 
ticle that emerges from the coagulation either starts its 
waiting time at the moment of the reaction ("kill") or 
continues evolving according to its prior setting, which 
is unchanged by the reaction ("no- kill"). The simula- 
tion then continues until it reaches the time t of interest 
or the concentration c(t) of interest. The experiment is 
then repeated over and over again for an ensemble of such 
chains. 

As mentioned earlier, there are issues related to 
the choice of the waiting time distribution, specifically 
whether to use the Pareto form Eq. ([4"4")l or the Mittag- 
Leffler form Eq. g5]). [UGJ]. For a < 1 they both give 
the same asymptotic results. At short times for a < 1 
and for all times when a = 1 the Mittag-Leffler distri- 
bution is the appropriate one to use. Most of the results 
presented below are asymptotic and for a < 1, and we 
mostly use the Pareto distribution, having ascertained 
to our satisfaction that both lead to the same outcome. 
For short-time results we use the Mittag-Leffler distribu- 
tion. Note that when looking at short-time behavior we 
need to specify the initial distribution of particles over 
the line. We have consistently chosen a random Poisson 
initial distribution. 



A. Results 

We will now test our theory against simulation results, 
and anticipate that the agreement is gratifyingly good. 
Figure [T] contains a variety of results for the concentra- 
tion c(t) of surviving A particles as a function of time 
for a = 0.5 and an initial concentration c(0) = 1. The 
initial concentration is sufficiently high for there to be 
essentially no transient behavior before the asymptotic 
power law dependence of c(i) takes over, as evidenced by 
the straight lines. First, note the comparison between 
the simulation (sim) results with a Pareto-like waiting 
time distribution (squares) and those of a Mittag-Leffler 
form (circles). Both are essentially linear, as expected. 
The slopes are the same, but the Pareto results are a 
little higher. A best fit leads to c PiSim ~ 0.751 i-°- 252 
and cm, sim ~ 0.541 t~ a25 °. The open symbols are "kill" 
simulations (and the fits just given are for this case), 
while the solid symbols are for "no-kill." The "no-kill" 
rule leads to consistently higher concentrations. This 
makes sense since the survivor at each reaction event is 
the A that arrived first at the reaction site, and it re- 
mains there longer (due to aging) than would the other 
reaction partner. This leads to its longer survival. 

While these results and comparisons are interesting, 
our most important task is to compare these results with 
those of our theory, which is shown by the upper solid 
line for the Pareto case and the lower solid line for the 
Mittag-Leffler case. The slopes are in excellent agree- 
ment with those of the simulations. The coefficients 
fall precisely between the "kill" and "no- kill" results. 
Specifically, we find that cp^theory ~ 0.819 i -0 - 25 and 
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FIG. 1: Reactant concentration vs time on a lattice of 10 4 
sites with periodic boundary conditions. c(0) = 1, a = 
0.5. Symbols give results of numerical simulations with kill 
(open symbols) and no-kill (black symbols) protocols. Cir- 
cles: Mittag-Leffler-based waiting time distribution. Squares: 
Pareto waiting time distribution. Solid lines: our theory, 
Eq. (|40)l . with K alpha . P (upper line) and with K a:M L (lower 
line). Broken lines: approximations c(t) ~ 1/Sp(t) (upper) 
and c(t) ~ l/SMh{t) (lower). For a definition of symbols and 
for quantitative fits to the various lines see text. In addition, 
in all the figures of this paper the error bars are smaller than 
the symbol sizes. 



CML.theory ~ 0.615 i _a25 . One might be tempted to 
conjecture that the theoretical model in some sense lies 
between the "kill" and "no-kill" scenarios. For example, 
perhaps the model is most germane if "kill" or "no-kill" 
are randomly selected according to some appropriate dis- 
tribution. For now, this remains a matter of conjecture. 

Finally, one additional result is shown in Fig. [TJ 
namely, the validity of the relation between the concen- 
tration of surviving reactant and the distinct number of 
sites visited by a random walker, c(i) ~ 1/S(t) @, [IB]- 
The dashed lines show 1/S(t), with the upper dashed line 
corresponding to the Pareto case, 1/Sp(t) ~ 0.853 <~ - 25 , 
and the lower dashed line associated with the Mittag- 
Leffler choice, \/S M L{t) ~ 0.641 t~°- 25 . 

Figure [5] again shows the concentration of surviving 
reactant as a function of time for a — 0.5 on a lattice 
of 10 sites, but now we explore whether our full theory, 
Eq. (|40p. in fact also captures the transient behavior. We 
explore this behavior by starting with a lower initial con- 
centration, c(0) = 0.1. Clearly, the theory again works 
extremely well for all times. Here the simulations are 
carried out in the "kill" scenario only, and we have pre- 
sented both the Pareto (squares) and Mittag-Leffler (cir- 
cles) results. The solid curves are the theoretical results 
for the Pareto (upper) and Mittag-Leffler (lower) cases. 
The dashed lines are the asymptotic results Eq. (13~1]) for 
the two cases. 

We stress that the asymptotic exponent of time as 
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FIG. 2: Reactant concentration vs time on a lattice of 10 4 
sites with periodic boundary conditions. c(0) =0.1, a = 0.5, 
showing transient behavior. Symbols give results of Numeri- 
cal simulations with "kill" protocol for Pareto (squares) and 
Mittag-Lefner-based (circles) waiting time distribution. Solid 
lines: our theory, Eq. (|40[) , with K a ,p (upper line) and K a ,ML 
(lower line). Broken lines: asymptotic approximation for c(t) 
as given in Eq. (I31|l for both cases. 

predicted by all theories (including our earlier theory 
that suffers from other difficulties) is correct and agrees 
with simulation results, which also agree with each other. 
The different theories and simulations (Pareto vs Mittag- 
Leffier, "kill" vs "no-kill," our earlier theory vs our cur- 
rent theory, distinct number of sites visited predictions) 
lead to different (but not wildly different) prefactors. A 
more stringent test is provided by the interparticle dis- 
tribution function, for which we now present a series of 
figures. 

In the ( "normal" ) reaction-diffusion problem the inter- 
particle distribution function develops a growing gap at 
small distances. This gap arises as spatially close pairs 
react and are not replenished because diffusion is slow in 
one dimension. The gap explains the "anomalous" decay 
law c(i) ~ i -1 / 2 , called anomalous because it differs from 
the law of mass action behavior c(t) ~ predicted for a 
Poissonian distribution of well-mixed reactants. Figure[3] 
shows these results as obtained from theory and simula- 
tions. They agree extremely well, which is a confirmation 
that we are in the asymptotic regime at time t = 6000. 
The exact theoretical expressions are well known [l2|, [l3[ , 

p(z) = ^ze-** 2 / 4 (46) 

and 

E(z) = erfc(V^z/2). (47) 

This figure serves as a basis of comparison for subsequent 
subdiffusive results. 

Figure Q] shows the same three panels for the subdif- 
fusive case a — 0.7, but it is necessary to go to longer 
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FIG. 3: Asymptotic interparticle distribution function and 
empty interval function vs the scaled variable z = c(t)x for 
the classic reaction-diffusion problem. The simulation results 
(symbols) are obtained with an exponential waiting time dis- 
tribution on a lattice of 40,000 sites at time t = 6000. The 
solid curves from the traditional exact theories are given in 
Eqs. (|46[1 and (|47[1 . The upper panel shows the interparticle 
distribution function and the inset shows it on a logarithmic 
scale that exhibits the simulation scatter at extremely low 
densities. The lower panel shows the empty interval function. 



times, t = 10 6 , to arrive at the asymptotic behavior. The 
theoretical curves are obtained from Eqs. and ([50]) 
and the simulations are carried out using a Pareto wait- 
ing time distribution with a "kill" rule. The agreement 
between theory and simulations is still excellent, with 
very small differences apparent near the maximum of the 
interparticle distance distribution. We have ascertained 
that we are in the asymptotic regime, so these differences 
are an indication of the approximate nature of the the- 
ory rather than of uncertainties pre-asymptotic transient 
effects. Figure [5] shows the three panels once again, 
now for the extremely subdiffusive case a = 0.2. The 
particles spend a great deal of time simply waiting be- 
tween jumps, and so it is necessary to go to much longer 
times, here t = 10 19 , to reach asymptotic behavior. The 
differences between theory and simulation are still small 
but certainly more noticeable. It is of course not clear 
whether the differences arise from the fact that the dis- 
tribution of particles is affected by the reaction, or from 
the assumption of linearity of the evolution operator of 
the empty interval function, or both. In any case, it does 
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FIG. 4: Asymptotic interparticle distribution function and 
empty interval function as a function of the scaled variable 
z = c(t)x for a — 0.7. The simulation results (symbols) 
are obtained with a Pareto waiting time distribution and a 
"kill" rule on a lattice of 40,000 sites at time t = 10 6 . The 
solid theoretical curves are obtained from Eqs. (|33p and (|30l) . 
The upper panel shows the interparticle distribution function 
and the inset shows it on a logarithmic scale that exhibits 
the simulation scatter at extremely low densities. The lower 
panel shows the empty interval function. 



FIG. 5: Asymptotic interparticle distribution function and 
empty interval function as a function of the scaled variable 
z = c(t)x for a — 0.2. The simulation results (symbols) 
are obtained with a Pareto waiting time distribution and a 
"kill" rule on a lattice of 40,000 sites at time t = 10 19 . The 
solid theoretical curves are obtained from Eqs. (|33|l and (|30[) . 
The upper panel shows the interparticle distribution function 
and the inset shows it on a logarithmic scale that exhibits 
the simulation scatter at extremely low densities. The lower 
panel shows the empty interval function. 



not seem an exaggeration to assert that the theory is very 
good. It is also clear that a gap in the interparticle dis- 
tance distribution develops no matter how subdiffusivc 
the system, contrary to our earlier thinking [5[. 

Two further issues are addressed in the following fig- 
ures. In most of our simulation results we have used the 
"kill" rule whereby the walker that arrives at an already 
occupied site eliminates the particle that was there, and 
we have stated that agreement with our theory is bet- 
ter with this rule than with the "no-kill" rule whereby 
the newcomer is eliminated. Figure [5] shows three panels 
where this is illustrated for a = 0.2, 0.5, and 0.8 (the 
rule choice does not matter when a = 1) . We see that as 
a decreases the differences in the simulation results for 
these two cases increase, and that the theory is closer to 
the "kill" results. Note that in these figures instead of a 
final simulation time we report a final simulation concen- 
tration to take advantage of the fact that z depends on 
time only through the concentration and hence we avoid 
additional uncertanities. These differences point to the 
approximation inherent in the hypothesis that the reac- 



tion does not alter the spatial distribution of reactants. 
It does, although not by very much, especially if a is not 
extremely small. 

Finally, we briefly examine the approach of the inter- 
particle distribution to its asymptotic behavior by plot- 
ting p(x,t)/c(t) as a function of xc(t), the scaling vari- 
able. Figure [7] shows the time progression of the distribu- 
tion for a = 0.5 until its arrival at asymptotic behavior 
in the lower panel. It is preferable to use the Mittag- 
Leffler waiting time distribution for this progression as 
we have done in these simulations because of its advan- 
tages at early times compared to the Pareto distribution. 
The figure illustrates that the growth of the interparticle 
gap at early times is faster than at later times when it 
is determined by the scaling form, but that the system 
settles into its asymptotic form rather quickly. 
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FIG. 6: Interparticle distribution function p(z). Symbols: 
simulation results with a Pareto waiting time distribution. 
Upside down triangles (red in color): "kill" rule. Right side 
up triangles (blue in color): "no-kill" rule. Lattice size L — 
10000, c(0) = l,c{tfi„ai) = 0.0025. Upper panel: a = 0.2. 
Middle panel: a = 0.5. Lower panel: a = 0.8. The solid 
curves are calculated from Eq. (|33p . 



VI. CONCLUSIONS 

We have presented a new theory for the coagulation 
reaction A + A — ► A on a lattice where the A's perform a 
subdiffusivc continuous time random walk characterized 
by the subdiffusive exponent a. Our theory relies on the 
connection between continuous time random walks and 
fractional diffusion equations, and is based on two as- 
sumptions. The first hypothesis is that the evolution of 
the probability density for an empty interval of length x 
at time t in the presence of reactions is the same as the 
probability density that, in the absence of reactions, the 
distance between two subdiffusive particles that start at 
the same location at t — is x at time t. The two prob- 
ability densities are not equal because they obey differ- 
ent initial and boundary conditions. Only the evolution 
equations are assumed to be the same. The second hy- 
pothesis is that the (unknown) common equation of evo- 
lution for these two probability densities is linear. We 
have no analytic way of checking the validity of these 
assumptions, but the results of the theory are extremely 
close to those of numerical simulations in all cases tested. 




xc(t) 



FIG. 7: Interparticle distribution function for a = 0.5. Sym- 
bols: simulation results with a Mittag-LefHer waiting time 
distribution with C(0) = 0.1. Upper panel: t = 10 2 . Middle 
panel: t = 10 5 . Lower panel: t = 10 s . The solid curves are 
asymptotic, calculated from Eq. (|33[) . 

Some of the small differences are more likely ascribed to 
the second hypothesis, the assumption of linearity of the 
evolution operator, than the first, as explained in the 
earlier analysis. 

We explicitly calculated the global reactant concentra- 
tions as a function of time. This is not a very discrimi- 
nating measure of the quality of different theories; all the 
theories and implementations tested in this work lead to 
the same correct decay exponent of time, c(t) ~ t~ a / 2 , 
and only the prefactors vary, albeit not dramatically. The 
well known association of c(t) with the inverse of the 
distinct number of sites visited by a subdiffusive walker 
also gives the correct exponent and a prefactor within 
the range of our theory and simulation results. We also 
tested the early time transient behavior of c(t) predicted 
by our theory and find that agreement with simulation 
results is also excellent. 

A more stringent test of the theory is the interparticle 
distribution function, that is, the more "local" function 
p(x,t). Again, we find our comparisons to show that the 
theory captures the simulation results very well, even for 
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extremely subdiffusive systems. Most of the results that 
we present are asymptotic, but we also tested the validity 
of our theory in the transient regime, and found again 
that it works equally well. 

An issue that we discussed at some length concerns a 
choice that must be made at the microscopic (continuous 
time random walk) level used in the simulations but not 
at level of the mesoscopic fractional diffusion equation 
used in the theory. The choice concerns the particle that 
vanishes when two A particles land at the same location, 
the one that was there already or the one that just ar- 
rived. The results of simulations show the differences to 
be small but discernible in one dimension, which is evi- 
dence of aging [10(; we find that the agreement with our 
theory is slightly better with the first choice and per- 
haps would be even better with a combination of choices 
that allows for both possibilities with some probability 
distribution. We also discussed the choice of the wait- 
ing time distribution, an issue that becomes particularly 
important at early times and when a approaches unity. 

Finally, we point to our work in higher dimensions [6| . 
In that work we considered the same reaction as well as 
the annihilation reaction but on a lattice with traps of 
random depths. The depths of the traps are associated 
with mean sojourn times that are finite and distributed 
according to a power law. Using a mean field formula- 
tion, we approximated this system by one that has iden- 
tical fat-tailed waiting time distributions for each site, 
and show that in dimension d = 3 the mean field model 
already captures the reaction dynamics and distribution 
functions of the original random trap system. In the 



mean field model, the "kill" and "no-kill" scenarios give 
the same results. A particularly important result in that 
work is that the waiting time distribution is unaffected 
by the reaction. In our one-dimensional model we start 
with the translationally invariant system, but also pos- 
tulate, and show to lead to excellent results, an equiv- 
alence between evolution operators in the systems with 
and without reactions. We also point to the fact that 
the distribution develops a growing gap between parti- 
cles @. The growth of the interparticle gap at early 
times is faster than at later times when it is determined 
by a scaling form, but in any case there is a gap. This 
"anomalous" distribution leads to the "anomalous" de- 
cay kinetics of the concentrarion, as in the case of normal 
diffusion in one dimension. 

The analysis presented here can be adapted to the 
A + A — > reaction, but the quantities to be calculated 
there are different because gap lengths don't change con- 
tinuously. One concentrates instead on the parity of the 
number of particles in a given interval [H, H[ , from which 
the desired quantities can in turn be derived. This will 
be a subject for future work. 
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